Un peu de contexte

Cette pandémie, dans laquelle un virus totalement nouveau se propage au sein de la population mondiale présente de nombreux inconnus et les incertitudes sont très grandes.

Dès le début de l’épidémie, une des principales questions péoccupantes à propos de cette maladie émergente concerne son taux de mortalité.

Cette information est importante car elle permet de nous faire rapidement une idée des conséquences les plus dramatiques que cette maladie peut avoir lorsqu’elle passe de l’épidémie à la pandémie.

Très rapidement après le début de l’épidémie, nous avons commencé à recevoir, notemment en provenance de Chine, des données concernant le nombre de malades identifiés (testés) et le nombre de victimes du nouveau virus. La source actuelle qui fait consensus concernant la collecte l’oganisation et la synthèse de ces informations se trouve ici : Johns Hopkins University Dashboard

La fameuse formule

Vue comme ça, cette formule semble peu parlante (sachez également qu’elle peut s’écrire sous différentes formes). Nous n’allons pas la décrypter en détail ici, nous nous contenterons simplement de noter quelques éléments nécessaires à nous éclairer sur les raisons de l’utilité de cette formule :

• P(A) signifie probabilité de A
• | signifie sachant
• P(A|B) signifie donc probabilité de A sachant B

Ces quelques détails fournis nous pouvons constater que la formule permet globalement de connaître la probabilité de A sachant B sur la base de la probabilité de A, de la probabilité de B et de la probabilité de B sachant A. C’est la raison pour laquelle Laplace parlait de probabilité inverses ;)

Si nous remplaçons A par taux de létalité et B par données nous voyons que notre formule nous permet de connaître la probabilité du taux de létalité sachant nos données sur la base de la probabilité des données sachant le taux de létalité.

Cette version de la formule nous indique simplement que la probabilité d’un taux de létalité sachant nos données est proportionnel au produit de la probabilité des données sachant ce taux et de la probabilité du taux.

Cette formule indique donc simplement qu’en combinant notre a priori sur le taux de létalité avec la probabilité que nos données puissent être générées dans un monde où un tel taux est vrai nous permet de calculer la probabilité de ce taux sachant nos données, c’est à dire notre a posteriori (pour plus de commodité nous parlerons aussi de prior pour désigner l’a priori et de posterior pour désigner l’a posteriori).

Implémenter un modèle génératif

Certes nous avons une description très générale du principe du modèle génératif, mais en pratique comment faire ??

Prenons un exemple, illustrons le et implémentons le pour nous en faire une meilleure idée.

Prenons le cas d’un foyer épidémique du MERS-CoV dont le taux de létalité sous-jacent serait de 30%. Imaginons maintenant que nous voulions simuler une situation où un épidémiologue se rend dans ce foyer épidémique et identifie qu’il comporte 160 malades. Quel nombre de décès pourrions nous envisager ? Bien entendu nous pourrions simplement calculer 30% de 160 malades (i.e. 48 victimes) et nous satifaire du résultat.

Mais nous savons que 30% est une sorte de valeur de la proportion de morts autour de laquelle notre donnée simulée doit se situer mais, l’aléa fait qu’il ne s’agit pas d’une valeur exacte et absolue en tout temps. Nous allons donc utiliser une distribution de probabilité dont le paramètre sera 30% et qui nous permettra de fabriquer une donnée du nombre de victimes pour 160 malades virtuels.

Notre modèle génératif va donc implémenter cette distribution de probabilité. Dans notre cas la sortie du modèle génératif doit nous indiquer pour chacun des 160 malades virtuels s’ils sont des victimes ou non. La distribution de probabilité capable de remplir ces critères est la distribution binomiale.

Le schéma suivant illustre le process que nous allons mettre en place pour générer une donnée simulée

En pratique le code R pour ce cas particulier de modèle génératif est le suivant :

victimes_virtuelles <- rbinom(1, size = 160, prob = 0.30)
victimes_virtuelles

## [1] 37

Si vous exécutez plusieurs fois de suite cette ligne de code vous constaterez qu’a chaque nouvelle exécution le nombre de victimes virtuelles change et gravite plus ou moins autour de la valeur 48.

Nous avons donc été capable de simuler un nombre de victimes virtuelles sur la base d’un taux de létalité donné.

Remonter à la valeur de paramètre

Notre problème c’est que ce qui nous intéresse c’est l’exact inverse !! Nous avons des données et nous cherchons à remonter au paramètre de taux de létalité, en d’autres termes plus techniques nous cherchons à faire une inférence bayésienne.

Cette petite vidéo de la chaîne Hygiène mentale expose ce même problème avec des lancés de dés. Nous savons qu’un dé a été lancé et qu’il a produit un 6, ce que nous cherchons à savoir : quel type de dé a été lancé ? Est-ce un dé à 4, 6, 8, 12 ou même 20 faces ? Comme pour notre taux de létalité, il faut procéder à une inférence bayésienne pour déterminer quel dé a le plus probablement été lancé.

Nous allons voir que malgré cette difficulté apparente, nous sommes sur le bon chemin pour passer du Prior au Posterior.

Implémentation du prior

En pratique on procède d’abord au tirage des 10000000 de valeurs de taux

Nous avons défini un Prior non informé, pour lequel chaque taux de létalité peut être équiprobablement situé entre 0 et 100%, nous pouvons l’implémenter comme suit :

nb_tirages <- 10000000
prior <- runif(n=nb_tirages, min=0, max=1)

Dans cette ligne de code, nous créons une liste appelée prior dans laquelle nous plaçons les résultats des 10000000 de tirages de valeurs comprises entre 0 et 1 (1 équivalant à 100%).

Nous pouvons l’illustrer graphiquement ainsi :

hist(prior, xlim=c(0,1), ylim=c(0,600000), 
     main="A Priori", xlab="Taux de létalité", 
     ylab="Fréquence", col="grey")

Nous constatons bien sur cet histogramme que chaque valeur de taux de létalité est tirée en quantité quasi-égale.

Implémentation du modèle génératif

Ensuite on fournit toutes ces valeurs de taux au modèle génératif.

Implémentons d’abord celui-ci

modele_generatif <- function(taux, nb_cas) {
  nb_morts <- rbinom(1, size = nb_cas, prob = taux)
  nb_morts
}

Cette fonction très simple prend en entrée le taux de létalité et le nombre de cas pour lequel on doit générer des données de nombre de morts. En sortie, elle retourne le nombre de morts simulés.

Simulation des données

Nous pouvons donc maintenant fournir tous nos tirages de taux au modèle génératif afin qu’il produise les données simulées :

data_cas = 271
morts_simu <- rep(NA, nb_tirages)
for(i in 1:nb_tirages) {
  morts_simu[i] <- modele_generatif(taux = prior[i], nb_cas = data_cas)
}

Ce petit morceau de code produit la liste des 1000000 données simulées.

Nous allons donc pouvoir filtrer ces données pour ne conserver que celles qui corresondent à notre donnée réelle : 1 décès.

posterior <- prior[morts_simu == 1]

Nous obtenons une liste de toutes les valeurs de taux qui parmi les 10000000 valeurs équiprobables tirées de notre distribution Prior ont abouti, lorsqu’elles ont été fournies au modèle génératif, à des données correpondantes à notre donnée de létalité pour le jour 1.

Nous devons vérifier que ce nombre de données simulées correspondant à nos donnés réelles soit suffisant.

length(posterior)

## [1] 36964

Cette liste n’est autre que notre distribution des taux a posteriori. Nous pouvons la représenter au moyen de l’histogramme suivant

hist(posterior,  breaks=40, axes=F,
     main="A posteriori", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))
axis(2)
axis(1, at=seq(0, 0.05, 0.001))

Nous pouvons à partir de cette distribution faire une estimation descriptive de l’intervalle de taux de létalité dans lequel se trouve 95% des taux de létalité qui ont été à l’origine de données simulées compatibles avec les notres et qui ont la plus grande fréquence d’apparition.

Techniquement cet intervalle traduit ce qui est appelé Intervalle de crédence à 95%

quantile(posterior, c(0.025, 0.5, 0.975))

##         2.5%          50%        97.5% 
## 0.0009101049 0.0061166739 0.0202749722

mean(posterior)

## [1] 0.007302888

On constate ainsi que 95% des taux qui ont produit des données compatibles se situent dans la fourchette de 0.08% et 2% avec un taux médian de 0.6% et un taux moyen de 0.7%

L’apport de nos données ont donc très clairement modifié notre vision du taux de létalité potentiel de la maladie. Nous sommes passé d’une situation où ce taux pouvait équiprobablement prendre n’importe quelle valeur entre 0 et 100% à une valeur possiblement comprise entre 0.08% et 2% en considérant nos données au jour 1.

La figure suivante représente notre passage du prior au posterior (attention les échelles des deux histogrammes sont très différentes ici)

Prudence

Cette prédiction est intéressante mais n’est-il pas un peu risqué de nous hasarder à des prédictions alors que nos estimations ne sont basées que sur la prise en compte des informations issue de 271 cas identifiés et 1 mort ?

Bien évidemment la réponse est oui, mais notons tout de même que l’approche est intéressante. et qu’elle nous permet de produire des anticipations avec peu de données et surtout, une certaine quantification de l’incertitude.

Avant de commencer à prendre en compte les données des jours suivant, attardons nous quelques temps sur l’effet des données que nous prenons en compte sur le posterior.

Village A

Reprenons le code exécuté précédemment en remplaçant les données du premier jour en Allemagne par les données du village A. Ici toutes les étapes de simulation des données sont les mêmes (morts_simu) seul le filtrage des données simulées est différent : là où nous ne retenions que les données simulées pour lesquelles le nombre de décès était 1 maintenant nous n’allons retenir que les données simulées pour lesquelles le nombre de décès est 135.

posterior_vill_A <- prior[morts_simu == 135]

length(posterior_vill_A)

## [1] 36743

hist(posterior_vill_A,  breaks=40, axes=T, xlim=c(0, 1), 
     main="A posteriori Village A", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.4))

quantile(posterior_vill_A, c(0.025, 0.5, 0.975))

##      2.5%       50%     97.5% 
## 0.4387548 0.4981299 0.5570957

mean(posterior_vill_A)

## [1] 0.498066

Comme Vous pouvez le constater, cette fois, le taux estimé dans le posterior est compris entre 43% et 55%. Nous confirmons donc bien que pour un tirage de nombre de malade identique, le nombre de décès influe correctement sur l’estimation du taux a posteriori

Pour rappel voici l’histogramme de la distribution du taux de létalité pour les données du premier jour en Allemagne

hist(posterior,  breaks=40, axes=T, xlim=c(0, 1), 
     main="A posteriori (prior uniforme)", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))

Voyons les données pour le village B

Village B

Cette fois les données à notre disposition sont bien moins nombreuses. Rerprenons le code pour observer l’effet de cette situation particulière.

Cette fois, il nous faut générer de nouvelles données simulées en procédant à un tirage de 6 malades pour chaque taux de la distribution prior.

data_cas = 6
morts_simu_B <- rep(NA, nb_tirages)
for(i in 1:nb_tirages) {
  morts_simu_B[i] <- modele_generatif(taux = prior[i], nb_cas = data_cas)
}

Nous filtrons les données simulée compatibles (3 morts)

posterior_vill_B <- prior[morts_simu_B == 3]

length(posterior_vill_B)

## [1] 1428956

Voyons l’histogramme de la distribution des taux a posteriori.

hist(posterior_vill_B,  breaks=50, axes=T, xlim=c(0, 1), 
     main="A posteriori Village B", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))

Comparativement au village A notre distribution des taux probables est toujours centrée autour de 50% toutefois, la largeur de la distribution est bien plus grande.

quantile(posterior_vill_B, c(0.025, 0.5, 0.975))

##      2.5%       50%     97.5% 
## 0.1844881 0.5002234 0.8162549

mean(posterior_vill_B)

## [1] 0.5001771

Cette fois le taux de létalité a posteriori se situe possiblement entre 18% et 81%. Nous interprétons cette augmentation de l’intervalle de crédence à 95% comme une mesure de l’augmentation de notre incertitude concernant le taux de létalité réel.

A titre d’ultime comparaison, voyons ce que nous obtenons comme distribution a posteriori avec les données du village C

Village C

data_cas = 10000
morts_simu_C <- rep(NA, nb_tirages)
for(i in 1:nb_tirages) {
  morts_simu_C[i] <- modele_generatif(taux = prior[i], nb_cas = data_cas)
}

Nous filtrons les données simulée compatibles (5000 morts)

posterior_vill_C <- prior[morts_simu_C == 5000]

length(posterior_vill_C)

## [1] 971

Histogramme de la distribution des taux a posteriori.

hist(posterior_vill_C,  breaks=25, axes=T, xlim=c(0, 1),
     main="A posteriori Village C", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))

quantile(posterior_vill_C, c(0.025, 0.5, 0.975))

##      2.5%       50%     97.5% 
## 0.4903229 0.5005129 0.5098661

mean(posterior_vill_C)

## [1] 0.5002361

Nous observons ici, que de données en très grande quantité, ont pour effet de réduire l’intervalle de crédence des taux les plus probables. En d’autres termes, plus nous intégrons de données à notre analyse, moins notre incertitude sera grande.

C’est l’un des grands points forts de l’analyse bayésienne de données, elle permet non seulement d’effectuer des estimations à partir de quantités de données variables même lorsqu’elles sont faibles et de surcroit, elles apportent une information sur l’incertitude associée à cette estimation.

Transformer le posterior en prior

Au moment où nous avons défini notre premier prior nous avons décidé qu’il s’étalerait entre 0 et 100% avec une équiprobabilité pour chaque taux. Mathématiquement il s’agit d’une distribution de probabilité dite uniforme.

Cette distribution est purement mathématique, de fait lors de l’implémentation du prior dans l’algorithme ABC, nous avons généré une liste de 10000000 de valeurs de taux en respectant cette loi mathématique. Notre véritable prior était donc une loi mathématique, que nous avons concrétisée au moyen de la création de cet énorme échantillon.

À la fin de l’exécution de notre algorithme nous avons obtenu un posterior constitué de la liste des taux qui avaient permis au modèle génératif de simuler des données identiques à nos données réelles. Ce que nous avons donc obtenu et illustré au travers d’histogrammes est un échantillon de valeurs de taux. Bien que pour des raisons pédagogiques nous ayons qualifié cet échantillon de posterior. En réalité il ne s’agit pas du véritable posterior. En effet, le véritable posterior est, comme le prior, une loi de probabilité. Un objet mathématique que nous avons approximé au travers de notre échantillon a posteriori.

Dans la mesure ou cet échantillon de taux posterior est suffisamment grand, nous pouvons considérer qu’il reflète la loi de probabilité que nous cherchions.

Pour réutiliser notre posterior en prior il va nous falloir trouver, à partir de cet echantillon posterior, quelle est cette loi de probabilité. Nous pouvons d’ores et déjà constater qu’il ne s’agit plus d’une loi uniforme. L’idéal serait de disposer d’une loi relativement flexible et il se trouve qu’il en existe une qui pourrait convenir : la distribution Beta.

La beta distribution comporte 2 paramètres α et β dont les valeurs permettent d’adapter la forme de la distribution.

Un autre avantage de cette distribution est qu’elle est bornée entre 0 et 1.

Regardons comment créer un échantillon d’1000000 de valeurs suivant cette loi de probabilité et faisons varier les paramètres α et β pour constater la flexibilité de cette loi.

α = 1 et β = 1

beta_sample <- rbeta(n = 1000000, shape1 = 1, shape2 = 1)

hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 1 et β = 1", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

Lorsqu’α = 1 et β = 1 la distribution Beta prend la forme d’une distribution uniforme

Essayons d’autres valeurs

• α = 1 et β = 5

beta_sample <- rbeta(n = 1000000, shape1 = 1, shape2 = 5)
hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 1 et β = 5", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

• α = 45 et β = 55

beta_sample <- rbeta(n = 1000000, shape1 = 45, shape2 = 55)
hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 45 et β = 55", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

• α = 7 et β = 3

beta_sample <- rbeta(n = 1000000, shape1 = 7, shape2 = 3)
hist(beta_sample, xlim=c(0, 1),
     main="Prior α = 7 et β = 3", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

• α = 50 et β = 50

beta_sample <- rbeta(n = 1000000, shape1 = 50, shape2 = 50)
hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 50 et β = 50", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

Nous avons donc à notre disposition une loi de probabilité dont les coefficients peuvent rendre compte de la distrubution des valeurs de taux posteriori que nous avons pu obtenir. Il nous faut maintenant reprendre notre échantillon de taux a prosteriori et déterminer les valeurs de coefficients α et β de sorte que nous puissions reprendre notre véritable posterior et le placer en tant que prior mis à jour.

Pour cela nous allons utiliser une fonction d’estimation des paramètres de la loi beta à notre distribution posterior et en récupérer les valeurs.

Cette fonction est fournie par le package EnvStats. La syntaxe est la suivante :

library(EnvStats)
beta_model <- ebeta(beta_sample, method = "mle")
coeff_alpha <- beta_model$parameters[[1]]
coeff_beta <- beta_model$parameters[[2]]

Dans cet exemple nous essayons de retrouver les coefficients de la dernière distribution beta que nous avons testés : α = 50 et β = 50

Nous obtenons ces coefficients

round(coeff_alpha)

## [1] 50

round(coeff_beta)

## [1] 50

Tout semble donc bien fonctionner.

Voyons quels sont les coefficients de la distribution beta qui correspondent à notre échantillon de taux posterior

post_beta_model <- ebeta(posterior, method = "mle")
coeff_alpha <- round(post_beta_model$parameters[[1]])
coeff_beta <- round(post_beta_model$parameters[[2]])
coeff_alpha

## [1] 2

coeff_beta

## [1] 270

Vérifions si la distribution de probabilité beta avec les coefficients trouvés correspond bien à celle que nous avons obtenus dans notre échantillon posterior

new_prior <- rbeta(n = 10000000, shape1 = coeff_alpha, shape2 = coeff_beta)
hist(new_prior, breaks = 60, axes=F, 
     main="Prior mis à jour", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(1,0,0,0.5))
axis(2)
axis(1, at=seq(0, 0.05, 0.001))

Tout semble donc en ordre !!

Nous avons bien déterminé notre nouveau prior !!

En procédant à cette vérification nous avons créé un nouvel échantillon de taux conformes à la loi beta de notre nouvel a priori.

ABC ne suffit pas, passons à JAGS

Nous pourrions recommencer comme nous l’avons fait jusqu’à présent, toutefois, si vous avez essayé de mettre en oeuvre les lignes de codes précédentes, vous vous êtes sans doute apperçu que son exécution prenait beaucoup de temps, voire même beaucoup trop longtemps si votre ordinateur manque un peu de ressources.

L’algorithme ABC n’est pas du tout efficace du point de vue computationnel, à partir de maintenant, nous allons employer des méthodes plus efficaces et donc plus rapides.

Nous ne nous attarderons pas ici sur le fonctionnement de ces algorithmes, sachez simplement que bien que fonctionnant différemment, ils suivent exactement les mêmes principes que ce que nous avons fait jusqu’à présent.

R n’est pas toujours efficace pour effectuer les simulations dont nous avons besoin, toutefois, il est très pratique pour manipuler et visualiser nos données. Nous allons donc utiliser un autre langage plus adapté que nous allons appeler directement via R.

Il existe de nombreux langages adaptés aux analyses bayésiennes de données. On les appelle communément les langages de programmation probabilistes. Parmi les plus connus on trouve : STAN, JAGS, WEBPPL, PYMC et bien d’autres.

Nous allons utiliser JAGS mais sachez que n’importe quel autre langage conviendrait tout aussi bien.

install.packages(rjags)

Première mise à jour

Dans un premier temps il nous faut implémenter le modèle en lui indiquant le prior et la loi qui dirige le modèle génératif. Nous avons établi notre distribution prior dans les sections précédentes et nous savons que la distribution de notre modèle génératif est une distribution binomiale.

L’implémentation impose trois étapes :

• La définition du modèle
• La compilation du modèle
• La simulation des données

library(rjags)

#DEFINE the model
covid_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

# Real Data
cases <- 802
deaths <- 2

#COMPILE the model
covid_jags <- jags.model(textConnection(covid_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
covid_sim <- coda.samples(model=covid_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

Visualisation du posterior

Nous pouvons donc aller droit au but et visualiser la distribution des taux Posterior issue de notre analyse par intégration des données du jour 2.

Pour une meilleure lecture nous avons ajouté la distribution de probabilité du prior en rouge puis superposé la distribution du posterior en vert.

# Create a df with covid_sim chains and another with the prior
covid_chains<- data.frame(covid_sim[[1]], iter = 1:10000)
df_prior_2 <- data.frame(new_prior)
colnames(df_prior_2)=c("p")

#PLOT the simulated posterior
ggplot(covid_chains, aes(x = p)) + 
  geom_density(data=df_prior_2, fill="red", alpha=0.4) +
  geom_density(fill="green", alpha=0.4) +
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Prior vs posterior jour 2')+
  theme(panel.background = element_blank())

Commme vous pouvez le constater, ces nouvelles données ont modifié la distribution de probabilité du taux de létalité.

Voyons quel est le nouvel intervalle de crédence à 95%

quantile(covid_chains$p, c(0.025, 0.5, 0.975))

##        2.5%         50%       97.5% 
## 0.001007512 0.003415019 0.008163162

Nous pensons donc désormais que le taux de létalité par cas est probablement compris entre 0.09% et 0.8%.

Nous avons bien procédé à une mise à jour bayésienne de notre a priori puis déterminé le nouveau posterior en incorporant les nouvelle données.

L’effet du prior

Nous avons précédemment consacré un temps à l’étude de l’effet des données sur le postérior lors de l’inférence bayésienne. Mais nous n’avons pas encore pris le temps d’étudier celui du prior. Or maintenant que nous avons mis à jour notre postérior sur la base d’un nouveau prior, il convient de mieux comprendre quel est l’effet de ce prior.

Nous allons donc reprendre notre exemple fictif avec un village D dans lequel les données sont de 50 cas et 25 morts (taux de létalité 50%). Cette fois au lieu de faire varier ces données, nous allons faire varier le prior.

Pour chaque graphique nous affichons en rouge le prior utilisé, en vert le posterior obtenu et en gris les précédents posteriors obtenus avec les précédents priors testés.

Avec un prior non informatif

• α = 1 et β = 1

coeff_alpha <- 1
coeff_beta <- 1
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)

#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(fill="green", alpha=0.4)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior non informatif')+
  theme(panel.background = element_blank())

old_post_vill_D <- vill_D_chains

Réitérons cette analyse en prenant différents priors :

• α = 1 et β = 5

coeff_alpha <- 1
coeff_beta <- 5
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")


#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(data=old_post_vill_D, fill="grey", alpha=0.4) +
  geom_density(fill="green", alpha=0.4)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior α = 1 et β = 5')+
  theme(panel.background = element_blank())

old_post_vill_D_2 <- vill_D_chains

• α = 45 et β = 55

coeff_alpha <- 45
coeff_beta <- 55
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")


#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(data=old_post_vill_D, fill="grey", alpha=0.2) +
  geom_density(data=old_post_vill_D_2, fill="grey", alpha=0.4) +
  geom_density(fill="green", alpha=0.4)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior α = 45 et β = 55')+
  theme(panel.background = element_blank())

old_post_vill_D_3 <- vill_D_chains

• α = 7 et β = 3

coeff_alpha <- 7
coeff_beta <- 3
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")


#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(data=old_post_vill_D, fill="grey", alpha=0.2) +
  geom_density(data=old_post_vill_D_2, fill="grey", alpha=0.4) +
  geom_density(data=old_post_vill_D_3, fill="grey", alpha=0.6) +
  geom_density(fill="green", alpha=0.8)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior α = 7 et β = 3')+
  theme(panel.background = element_blank())

Nous constatons ici que les mêmes données ne nous mènent pas au même posterior suivant le prior que nous employons. En effet, le posterior que nous obtenons est une forme de compromis entre le taux de létalité des données et le prior.

Comme pour la variation de la quantité des données, plus le prior est certain (distribution peu étalée), plus son effet sera fort. C’est à dire que le posterior obtenu sera proche de celui-ci. Ainsi dans une situation où le prior présente une très faible incertitude et les données sont présentes en faible quantité, le posterior sera très proche du prior. À l’opposé, si le prior est très incertain (par ex. prior non informatif) et les données en grande quantité, alors le posterior sera plus proche de la tendance centrale des données.

Poursuivons maintenant notre analyse.

Algorithme itératif

Pour cela nous allons mettre en place une boucle de répétition qui reprend automatiquement toutes les étapes que nous avons vues jusqu’à présent pour chaque ligne du tableau de données.

Pour mieux visualiser le résultat de notre algorithme, nous allons produire pour chaque jour un graphique de la distribution a posteriori.

Dans un premier temps, repartons du début en définissant les coefficients alpha et beta de la distribution a priori à la valeur 1 pour que nous partions d’un prior non informatif. Nous en profitons au passage pour indiquer sur quelle table de données nous souhaitons itérer.

Nous fabriquons aussi une table de données dans laquelle nous stockerons les intervalles de crédences et médiane de chaque distribution posterior.

## Empty df for fatality rate update tracing
rate_df = data.frame(day=integer(), Q25=double(), Q50=double(), Q97=double())

## Bayesian model preparation
# select the data table used to feed the process
data_table <- german_data

# Set initial Beta distribution coeffs correponding to a uniform distribution
coeff_alpha = 1
coeff_beta = 1

Puis, implémentons la boucle en charge du calcul itératif du modèle.

Voici une description de la liste des étapes de l’algorithme

1. Récupérer pour chaque ligne de la table de données parcourue, le nombre de cas identifiés et le nombre de décès enregistrés
2. Définir le modèle pour rjags
3. Compiler le modèle
4. Simuler 1000000 paramètres du posterior (Markov Chain Monte Carlo)
5. Stocker les simulations dans un data frame
6. Estimer les nouveaux paramètres α et β de la beta distribution pour la prochaine itération
7. Calculer et ajouter l’intervale de crédence et la médiane du posterior dans une table
8. Créer un graphique traduisant la nouvelle distribution du posterior

À la fin de la boucle nous ajoutons une dernière étape de création du graphique d’évolution de l’intervale de crédence et de la médiane du posterior au fil des jours

## Bayesian updating 
# For loop that iterates on the german-data table and proceed to the modelling steps
# update the beta distribution to set the daily new priors 
# and plot everyday death rate parameter probability distribution
for (row in 1:nrow(data_table)) {
  
  # 1. retrieve the cases and deaths daily quantities
  cases <- data_table[row, "cases"]
  deaths  <- data_table[row, "deaths"]
  
  # 2. Define the model
  covid_model <- "model{
              # Likelihood model for X
              X ~ dbin(p, n)

              # Prior model for p
              p ~ dbeta(a, b)
              }"
  
  # 3. Compile the model
  covid_jags <- jags.model(textConnection(covid_model), 
                           data=list(a=coeff_alpha, 
                                     b=coeff_beta, 
                                     X=deaths, 
                                     n=cases),
                           inits=list(.RNG.name="base::Wichmann-Hill",
                                      .RNG.seed=100))
  
  # 4. Simulate the posterior
  covid_sim <- coda.samples(model=covid_jags, 
                            variable.names=c("p"), 
                            n.iter= 1000000)
  
  # 5. Create a df with covid_sim chain
  covid_chain<- data.frame(covid_sim[[1]], iter = 1:1000000)

  
  # 6. Approximate the beta parameters  
  # (EnvStats package function)
  coeffs <-  ebeta(covid_chain$p, method = "mle")
  coeff_alpha <- coeffs$parameters[[1]]
  coeff_beta <- coeffs$parameters[[2]]
  

  # 7. Append posterior distribution 2.5%, 50% & 97.5% quantile values to
  # to fatality rate estimation table
  day_number=row
  chain_quant <- quantile(covid_chain$p, c(0.025, 0.5, 0.975))
  Q2 <- chain_quant[[1]]
  Q5 <- chain_quant[[2]]
  Q9 <- chain_quant[[3]]
  quantile_row <- cbind(day_number, Q2, Q5, Q9)
  rate_df <- rbind(rate_df, quantile_row)
  
  # 8. PLOT the simulated posterior
  day <- paste("Jour ", row)
  cas <- paste(cases, "malades")
  morts <- paste(deaths, "morts")
  a <- paste("alpha : ", coeff_alpha)
  b <- paste("beta : ", coeff_beta)
  plot <- ggplot(covid_chain, aes(x = p)) + 
              geom_density(fill="green", alpha=0.4, position = "stack")+
              xlim(0, 0.05) +
              ylim(0, 6000) +
              xlab('Taux de létalité') +
              ylab('Densité') +
              ggtitle(paste0(day, ", ",cas, ", ",morts, ", ",a, ", ",b))
  print(plot)
  file_name <- paste0(path, day, ".png")
  ggsave(filename = file_name)
  }
# End of the loop


# 9. Plot the posterior evolution over daily iterations
ggplot(rate_df, aes(x=rate_df$day, y=rate_df$Q5))+
  geom_line(aes(x=rate_df$day, y=rate_df$Q9), color="grey")+
  geom_line(aes(x=rate_df$day, y=rate_df$Q2), color="grey")+
  geom_ribbon(data=rate_df, 
              aes(ymin=rate_df$Q2,ymax=rate_df$Q9), 
                  fill="lightgrey", alpha="0.7")+
  geom_line()+
  ylim(0, 0.036) +
  xlab('Jours') +
  ylab('Taux de létalité') +
  ggtitle('Évolution de la distribution posterior du taux de létalité')+
  theme(panel.background = element_blank())
ggsave(filename = paste0(path, 'rate_evolution.png'))

Graphiques de sorties

Voici une animation de quelques un des graphiques produits par l’exécution de l’algorithme.

Ainsi que le graphique d’évolution du taux de létalité estimé a posteriori.

Interprétation

Nous pouvons faire deux constats principaux sur la base de ces deux représentations grahiques :

• Dans les premiers jours, l’incertitude concernant le taux de létalité est relativement grande mais diminue au fil du temps.

• Dans les jours suivants, l’incertitude reste globalement d’amplitude similaire mais le taux varie dans le sens d’une augmentation régulière.

Le premier constat s’explique en grande partie par le fait que le nombre de cas et le nombre de morts restent faibles dans les premiers jours.

Le second est plus problématique car il semble que notre procédure ne parvient pas à converger vers un taux de létalité donné. L’une des explications pourrait être liée au fait qu’il y a un délai entre le diagnostic de la maladie et le décès des patients (élément que nous avons volontairement mis de côté au début de notre analyse). Ce délai implique que les chiffres obtenus au début de la phase de croissance de l’épidémie contiennent une proportion plus importante de nouveaux malades que de nouveaux morts. Tous les malades étant récents, le délai pour qu’une fraction d’entre eux décède ne s’est pas encore écoulé. L’inverse se produit au moment de la phase de décrue de l’épidémie. En effet, lors du reflux de l’épisode épidémique, le nombre de nouveaux malades commence par décroître mais n’est pas immédiatement suivi par le nombre de décès qui peut ,temporairement, poursuivre sa croissance ou se stabiliser. À la fin de la vague épidémique on ne constate aucun nouveaux cas, mais une poursuite des décès.

Notre parti pris d’ignorer ce délai nous a donc ici porté préjudice. Pour autant, nous pouvons facilement imaginer des situations dans lesquelles un tel délai n’existe pas et donc dans lesquelles notre stratégie de mise à jour serait valide, comme des situations où des individus seraient amenés à prendre une décision oui / non par exemple.